home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
FishMarket 1.0
/
FishMarket v1.0.iso
/
fishies
/
376-400
/
disk_386
/
xlispstat
/
src2.lzh
/
XLisp-Stat
/
uni.c
< prev
next >
Wrap
C/C++ Source or Header
|
1990-10-02
|
3KB
|
166 lines
/* uni - uniform random number generator - Marsaglia's portable gen. */
/* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */
/* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */
/* You may give out copies of this software; for conditions see the */
/* file COPYING included with this distribution. */
#include <time.h>
#include "xlisp.h"
#include "osdef.h"
#ifdef ANSI
#include "xlproto.h"
#include "xlsproto.h"
#else
3include "xlfun.h"
#include "xlsfun.h"
#endif ANSI
#include "xlsvar.h"
#ifdef ANSI
double uni1(void);
void checkstate(LVAL);
int random_state_p(LVAL);
#else
double uni1();
void checkstate();
int random_state_p();
#endif ANSI
#define STATELENGTH 22
#define MDIG 32
#define M1 18
#define M2 19
#define I 20
#define J 21
/* seed: seed, m[1] - m[17], m1, m2, i, j */
LVAL newrandomstate()
{
LVAL state;
unsigned long seed;
long m1, m2, k0, k1, j0, j1;
int i;
xlsave1(state);
state = newvector(STATELENGTH);
/* set seed from clock */
seed = time((unsigned long *) NULL);
m1 = floor(pow(2.0, (MDIG - 2.0)) + .1)
+ floor(pow(2.0, (MDIG - 2.0)) + .1) - 1;
m2 = floor(pow(2.0, (MDIG / 2.0)) + .1);
setelement(state, 0, cvfixnum((FIXTYPE) seed));
setelement(state, M1, cvfixnum((FIXTYPE) m1));
setelement(state, M2, cvfixnum((FIXTYPE) m2));
setelement(state, I, cvfixnum((FIXTYPE) 5));
setelement(state, J, cvfixnum((FIXTYPE) 17));
if (seed % 2 == 0) seed--;
k0 = 9069 % m2;
k1 = 9069 / m2;
j0 = seed % m2;
j1 = seed / m2;
for (i = 1; i <= 17; i++) {
seed = j0 * k0;
j1 = ((seed / m2) + j0 * k1 + j1 * k0) % (m2 / 2);
j0 = seed % m2;
setelement(state, i, cvfixnum((FIXTYPE) j0 + m2 * j1));
}
xlpop();
return(state);
}
static double uni1()
{
LVAL state;
long k, i, j, m1;
state = getvalue(s_random_state);
checkstate(state);
m1 = getfixnum(getelement(state, M1));
i = getfixnum(getelement(state, I));
j = getfixnum(getelement(state, J));
k = getfixnum(getelement(state, i)) - getfixnum(getelement(state, j));
if (k < 0) k = k + m1;
setelement(state, j, cvfixnum((FIXTYPE) k));
i = i - 1;
if (i <= 0) i = 17;
j = j - 1;
if (j <= 0) j = 17;
setelement(state, I, cvfixnum((FIXTYPE) i));
setelement(state, J, cvfixnum((FIXTYPE) j));
return(((double) k) / m1);
}
double uni()
{
double x;
do { x = uni1(); } while (x <= 0.0 || x >= 1.0);
return (x);
}
int osrand(n)
unsigned n;
{
int k;
do {
k = n * uni();
} while (k < 0 || k >= n);
return(k);
}
LVAL xsmake_random_state()
{
LVAL arg = NIL;
if (moreargs()) arg = xlgetarg();
xllastarg();
if (arg == NIL) return(copyvector(getvalue(s_random_state)));
else if (arg == s_true) return(newrandomstate());
else {
checkstate(arg);
return(copyvector(arg));
}
}
LVAL xsrandom_state_p()
{
LVAL state;
state = xlgetarg();
return((random_state_p(state)) ? s_true : NIL);
}
static random_state_p(state)
LVAL state;
{
long i, j, m1;
if (! vectorp(state) || getsize(state) != STATELENGTH)
return(FALSE);
m1 = getfixnum(getelement(state, M1));
i = getfixnum(getelement(state, I));
j = getfixnum(getelement(state, J));
if (m1 <= 0 || i <= 0 || i > 17 || j <= 0 || j > 17)
return(FALSE);
return(TRUE);
}
static void checkstate(state)
LVAL state;
{
if (! random_state_p(state)) xlerror("bad random state", state);
}